setwd("/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison/NSAF/new analysis/SpC avgd across tech reps/")

OA.dat<-read.csv('GO slim xy plots for OA.csv', header=T)
#log transformation makes all numbers the same

#plot.nmds<-par(mfrow=c(3,4))

#SIMILAR EXPRESSION PROFILES BETWEEN TREATMENTS
#plot.small<-par(mfrow=c(1,3), mar=c(20,4,4,4))
x<-OA.dat$Org.Biogen.400[1:281]
y<-OA.dat$Org.Biogen.2800[1:281]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Cell Organization and Biogenesis', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:71], pts=c(13,236,69,144,219,61,57,248,227,220,249,215,132,90,244,169,225,176,51,33,188,160,79,77,80,129,9,62,99,35,71,104,122,136,117,5,182,19,16,37,133,4,28,218,91,225,119,86,226,53,162,264,247,217,127,85,280,7), col=c('green', rep('magenta',57)), pos=1)
#identify(x,y, labels=OA.dat$Org.Biogen[1:574], cex=0.5, n=2, col='magenta')
#few differences

x<-OA.dat$Cell.signaling.400[1:34]
y<-OA.dat$Cell.signaling.2800[1:34]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Cell-Cell Signaling', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:34], pts=c(1,18,6,33,32,9,27,2), col='magenta',pos=1)
#identify(x,y, labels=OA.dat$Cell.signaling[1:51], cex=0.5, n=5, col='magenta')
#small differences


x<-OA.dat$Protein.metabolism.400[1:386]
y<-OA.dat$Protein.metabolism.2800[1:386]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='Protein Metabolism', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:386], pts=c(27,56,85,50,53,140,87,62,326,322,237,137,375,25,161,19,47,306, 271,336,310,138,280,286,329,98,298,180,355,357,11,318,367,68,249,244,195,373,338,104,184,7,45,84,173,291,362,46,189,283,252,349,74,152,222,215,317,211,135,358,293,73,216,86,113,31,96,165,160,77,316,124), col='magenta',pos=1)
#small differences

#LARGE EXPRESSION DIFFERENCES
plot(OA.dat$Cell.adhesion.400, OA.dat$Cell.adhesion.2800, main='Cell Adhesion', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='magenta')
#differences
x<-OA.dat$Cell.adhesion.400[1:71]
y<-OA.dat$Cell.adhesion.2800[1:71]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Cell Adhesion', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1, col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:71], pts=c(46,39,2,18,45,6,4,23,35,31,47,48,71,68,41,70,49,38,37), col='magenta', pos=1)
#red dotted lines = 95% CI and blue line = 1:1 line
#highlight(y~x, data=df, lbls=OA.dat$Cell.adhesion[1:71], pts=c(46,39,2,18,45,6,4,23,35,31,47,48,71,68,41,70,49,38,37), col='black', cex=0.5)
#identify(x,y, labels=OA.dat$Cell.adhesion[1:109], cex=0.5, n=5)
#med-large differences

x<-OA.dat$RNA.metabolism.400[1:180]
y<-OA.dat$RNA.metabolism.2800[1:180]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='RNA Metabolism', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:180], pts=c(57,142,36,103,122,44,112,40,13,60,92,2,50,73,42,31,3,4,39,90,76,33,174,5,65,10,148,17,129,179,25,134,58,115,96,151,164,144,99,84,157,22,167), col='magenta', pos=1)
#big differences

x<-OA.dat$DNAmetabolism.400[1:32]
y<-OA.dat$DNAmetabolism.2800[1:32]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='DNA metabolism', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:32], pts=c(3,32,14,27,29,18,19,4), col=c('green',rep('magenta',7)), pos=1)
#identify(x,y, labels=OA.dat$DNA.metabolism[1:58], cex=0.5, n=7)
#large differences

#MODERATE DIFFERENCES
x<-OA.dat$Cycle.Prolif.400[1:103]
y<-OA.dat$Cycle.Prolif.2800[1:103]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Cell Cycle and Proliferation', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:103], pts=c(87,18,92,79,64,11,39,15,2,19,83,27,36,41,80,81,52,49,90,22), col='magenta', pos=1)
#medium differences

x<-OA.dat$Signal.trans.400[1:162]
y<-OA.dat$Signal.trans.2800[1:162]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='Signal Transduction', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:162], pts=c(44,31,142,1,33,92,60,125,34,78,46,8,17,24,49,5,4,26,2,68,137,159,79,88,73,39,22,6,148,149,102,136,141,15,114), col='magenta', pos=1)
#medium differences

x<-OA.dat$Stress.400[1:130]
y<-OA.dat$Stress.2800[1:130]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='Stress Response', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:130], pts=c(21,13,5,130,54,113,107,63,114,37,2,31,111,67,87,26,62,59,79,72,48,44,25,74,29,33,6,27,12,9), col='magenta', pos=1)
#medium differences

x<-OA.dat$Transport.400[1:260]
y<-OA.dat$Transport.2800[1:260]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df,  main='Transport', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:260], pts=c(72,130,58,252,66,61,186,208,240,192,48,214,161,126,168,59,189,8,116,6,24,73,37,57,132,78,41,181,158,136,80,83,127,7,4,97,9,20,101,249,14,113,146,3,172,32,47,100,121,16,22,182,42,25,171,241,92,36,87,244,103,196,223,159,156,213,235,245,43,230,152,178,74), col='magenta', pos=1)
#med differences

x<-OA.dat$Death.400[1:79]
y<-OA.dat$Death.2800[1:79]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='Death', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:79], pts=c(4,62,29,7,60,41,67,72,19,24,59,55,16,44,15,21,17,73), col='magenta', pos=1)
#medium differences

x<-OA.dat$Dev.processes.400[1:209]
y<-OA.dat$Dev.processes.2800[1:209]
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x,data=df, main='Developmental Processes', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')
lines(newx, preds[,3], lty='dashed', col='red')
lines(newx, preds[,2], lty='dashed', col='red')
highlight(y~x, data=df, lbls=OA.dat$Sig[1:209], pts=c(13,9,42,43,2,97,177,5,16,22,40,96,4,18,24,31,93,173,52,179,135,187,105,74,26,65,167,176,202,56,39,186,166,145,61,79,178,191,129,140,127,164,37,69,103,190,83), col=c('green', rep('magenta', 46)), pos=1)
#medium differences

#xy plots for all proteins
dat.xy<-read.csv('avg NSAF by treatment.csv', header=T, row.names=1)
library(NCStats)

#dat.xytra<-(dat.xy[,1:4]*10000)
#dat.xytra<-data.trans(dat.xytra, method='log', plot=F)
#sig<-dat.xy$Sig
#data.xytra<-cbind(dat.xytra, sig)

#plot(data.xytra$Avg.2800~data.xytra$Avg.400, xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
#highlight(data.xytra$Avg.2800~data.xytra$Avg.400, lbls=data.xytra$sig, pts=c(217,110,224,274,223,296,10,4,72,277,77,278,17,70,73,74,132,6,301,341,31,80,14,28,265,18,88,43,251,264,7,271,203,24,268,12,98,138,356,207,250,47,287,228,140,45,97,142,179,145,62,144,290, 68,107,362,108,364,365,366,367,368,369,370,1440,1441,1442,1443,1444,1445,1446,1447,1461,1535,1472,1491,1493,1615,1506,1605,1479,1485,1468,1514,1521,1602,1610,1473,1489,1571,1466,1476,1570,1464,1492,1471,1542,1451,1463,1522,1597,1482,1507,1568,1455,1486,1509,1454,1508,1604,1452,1480,1502,1518,1531,1458,1608,1528,1538), col='magenta', pos=1)
#legend(0, 3.25, col='magenta', legend='>5-fold difference', pch="*")

x<-dat.xy$Avg.2800
y<-dat.xy$Avg.2800MechS
plot(y~x, main='Mech Stress Response 2800', xlab='Expression at 2800 µatm', ylab='Expression at 2800 µatm + Mech Stress')
abline(0,1,col='blue')

x<-dat.xy$Avg.400
y<-dat.xy$Avg.400MechS
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Mech Stress Response 400', xlab='Expression at 400 µatm', ylab='Expression at 400 µatm + Mech Stress')
abline(0,1,col='blue')

x<-dat.xy$Avg.400
y<-dat.xy$Avg.2800
df<-data.frame(x=x, y=y)
mod<-lm(y~x, data=df)
newx<-seq(min(df$x), max(df$x), length.out=100)
preds<-predict(mod, newdata = data.frame(x=newx), interval='confidence')
plot(y~x, data=df, main='Ocean Acidification Response', xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm')
abline(0,1,col='blue')

setwd("/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison/NSAF/new analysis/SpC avgd across tech reps/go slim xy plots")
#data set only gives colors to proteins that are considered differentially expressed in manuscript: expressed in at least 2 oysters
dat.col<-read.csv('avg NSAF by treatment colors2.csv', header=T, row.names=1)
#dat.col<-read.csv('avg NSAF by treatment colors.csv', header=T, row.names=1)
col.tra<-(dat.col[,1:4]*10000)
col.tra<-data.trans(col.tra, method='log', plot=F)


col.oa<-as.character(dat.col$ColOA)
plot(Avg.2800~Avg.400, data=col.tra, col=col.oa, xlab='Expression at 400 µatm', ylab='Expression at 2800 µatm', pch=dat.col$pchOA, cex.lab=1.5)
legend(0, 3.25, col=c('forestgreen','magenta', 'blue', 'orangered3'), legend=c('Expressed in 1 treatment','>5-fold difference', '>2-fold difference', 'q-value < 0.1'), pch=c(rep(1,3),16))

col.400ms<-as.character(dat.col$Col400MS)
plot(Avg.400MechS~Avg.400, data=col.tra, col=col.400ms, xlab='Expression at 400 µatm', ylab='Expression at 400 µatm + Mech Stress', pch=dat.col$pch400MS, cex.lab=1.5)
legend(0, 3.25, col=c('forestgreen','magenta', 'blue', 'orangered3'), legend=c('Expressed in 1 treatment','>5-fold difference', '>2-fold difference', 'q-value < 0.1'), pch=c(rep(1,3),16))

col.2800ms<-as.character(dat.col$Col2800MS)
plot(Avg.2800MechS~Avg.2800, data=col.tra, col=col.2800ms, xlab='Expression at 2800 µatm', ylab='Expression at 2800 µatm + Mech Stress', pch=dat.col$pch2800MS, cex.lab=1.5)
legend(0, 3.25, col=c('forestgreen','magenta', 'blue', 'orangered3'), legend=c('Expressed in 1 treatment','>5-fold difference', '>2-fold difference', 'q-value < 0.1'), pch=c(rep(1,3),16))

col.ms<-as.character(dat.col$ColMS)
plot(Avg.2800MechS~Avg.400MechS, data=col.tra, col=col.ms, xlab='Expression at 400 µatm + Mech Stress', ylab='Expression at 2800 µatm + Mech Stress', cex.lab=1.5, pch=dat.col$pchMS)
legend(0, 3.25, col=c('forestgreen','magenta', 'blue', 'orangered3'), legend=c('Expressed in one treatment','>5-fold difference', '>2-fold difference', 'q-value < 0.1'), pch=c(rep(1,3),16))